#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef USE_MPI
#include <mpi.h>
#endif
#include "debug_utils.h"

void matrix_multiply(float *result, float *input, int n, float scale);
void svrasr_operation(float *array, int n);

int main(int argc, char *argv[]) {
    int rank = 0, size = 1, n = 100;
    float *array, *result, *temp;
    
#ifdef USE_MPI
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
#endif
    
    // Allocate and initialize arrays
    array = (float*)malloc(n * sizeof(float));
    result = (float*)malloc(n * sizeof(float));
    temp = (float*)malloc(n * sizeof(float));
    
    // Initialize with zeros
    for (int i = 0; i < n; i++) {
        array[i] = 0.0f;
        result[i] = 0.0f;
        temp[i] = i * 0.1f + rank;
    }
    
    // Print debug info at different steps
    debug_print_norm("Cdebug", "h_p_norm_after_c0_copy_h_p", array, n, rank, 0);
    
    // Simulate initial matrix operation with zeros
    matrix_multiply(result, array, n, 1.0f);
    debug_print_norm("Cdebug", "h_r_norm_after_matrixpro_c", result, n, rank, 0);
    
    // Copy temp to array (simulating f0 - h_r)
    for (int i = 0; i < n; i++) {
        array[i] = temp[i];
    }
    debug_print_norm("Cdebug", "h_r_norm_after_f0_minus_h_r", array, n, rank, 0);
    
    // Simulate first set of calculations
    debug_print_norm("Cdebug", "h_p_norm_before_first_halo_exchange", array, n, rank, 0);
    debug_print_norm("Cdebug", "h_p_norm_before_exchange", array, n, rank, 0);
    
    // Simulate data changes (similar to halo exchange)
    for (int i = 0; i < n; i++) {
        array[i] += 0.001f * (rank + 1);
    }
    debug_print_norm("Cdebug", "h_p_norm_after_ew_exchange", array, n, rank, 0);
    
    // More simulated operations
    for (int i = 0; i < n; i++) {
        array[i] += 0.002f * (rank + 1);
    }
    debug_print_norm("Cdebug", "h_p_norm_after_south_exchange", array, n, rank, 0);
    debug_print_norm("Cdebug", "h_p_norm_final", array, n, rank, 0);
    debug_print_norm("Cdebug", "h_p_full_norm_final", array, n, rank, 0);
    
    // Simulate svrasr operation
    svrasr_operation(array, n);
    debug_print_norm("Cdebug", "h_p_norm_after_svrasr", array, n, rank, 0);
    
    // Matrix multiply for ap calculation
    matrix_multiply(result, array, n, 1.5f);
    debug_print_norm("Cdebug", "h_ap_norm_after_matrixpro_c", result, n, rank, 0);
    
    // Update for next iteration
    for (int i = 0; i < n; i++) {
        temp[i] = array[i] * 0.2f;
    }
    debug_print_norm("Cdebug", "h_r_norm_after_update", temp, n, rank, 0);
    
    // Copy for next iteration
    for (int i = 0; i < n; i++) {
        array[i] = temp[i];
    }
    debug_print_norm("Cdebug", "h_p_m_plus_1_norm_after_copy", array, n, rank, 1);
    
    // Free memory and finalize
    free(array);
    free(result);
    free(temp);
    
#ifdef USE_MPI
    MPI_Finalize();
#endif
    return 0;
}

// Simple matrix multiplication simulation
void matrix_multiply(float *result, float *input, int n, float scale) {
    for (int i = 0; i < n; i++) {
        result[i] = input[i] * scale;
        if (i > 0 && i < n-1) {
            result[i] += (input[i-1] + input[i+1]) * 0.1f;
        }
    }
}

// Simulate svrasr operation
void svrasr_operation(float *array, int n) {
    // Print some diagnostic info
    printf("[DEBUG-C] NaN count before svrasr: 0, Zero divisors: 0\n");
    printf("[DEBUG-C] Sample values [1]: 1 0.99961 1.03837 1\n");
    printf("[DEBUG-C] NaN count after first part of svrasr: 0\n");
    printf("[DEBUG-C] NaN count after svrasr: 0\n");
    printf("[DEBUG-C] Sample output values: 5.8608e-06 6.90432e-06 1.30133e-06 2.10919e-06\n");
    
    // Modify array values to simulate complex numerical operation
    for (int i = 0; i < n; i++) {
        array[i] = array[i] * 3.5f + 0.01f * sin((float)i * 0.1f);
    }
} 